Baf53cre ENS neurons, SI data
Nb infection 5d, PBS(control) x4 INF(inflammation) x4
loading 8k cells for each,
demo called 20,985 cells
plus called 21,294 cells
cleaning-up and re-clustering and annotation
ref-mapping to NatNeur2021
mainly follow https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
GEX.seur <- readRDS("./GEX.seur.preAnno0717.rds")
GEX.seur
## An object of class Seurat
## 23258 features across 9423 samples within 1 assay
## Active assay: RNA (23258 features, 1500 variable features)
## 3 dimensional reductions calculated: pca, tsne, umap
color.FB <- ggsci::pal_d3("category20c")(20)[c(2,7,12,17,
1,6,11,16
)]
# Visualize QC metrics as a violin plot
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0.01, group.by = "FB.info", cols = color.FB)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0, group.by = "FB.info", cols = color.FB)
plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "FB.info", cols = color.FB)
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "FB.info", cols = color.FB)
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "FB.info", cols = color.FB)
plota + plotb + plotc
head(GEX.seur@meta.data)
## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCACAAGACGAC-1 Baf53Nb_Ileum 3257 1801 0.36843721 0.3991403
## AAACCCAGTACAGGTG-1 Baf53Nb_Ileum 5446 2500 0.09181050 0.2937936
## AAACCCAGTGGGCTCT-1 Baf53Nb_Ileum 1511 997 0.66181337 0.4632694
## AAACCCAGTTTGTTCT-1 Baf53Nb_Ileum 2855 1577 0.98073555 0.3152364
## AAACCCATCCTAGCCT-1 Baf53Nb_Ileum 2433 1451 0.08220304 0.3699137
## AAACCCATCGAAACAA-1 Baf53Nb_Ileum 3129 1656 0.12783637 0.4474273
## FB.info S.Score G2M.Score Phase RNA_snn_res.1.5
## AAACCCACAAGACGAC-1 Baf53Nb.IF4 0.011590405 -0.0004169865 S 3
## AAACCCAGTACAGGTG-1 Baf53Nb.CL3 -0.038183546 0.0022744719 G2M 0
## AAACCCAGTGGGCTCT-1 Baf53Nb.CL4 -0.024203070 0.0012414826 G2M 7
## AAACCCAGTTTGTTCT-1 Baf53Nb.IF1 -0.013980476 0.0039329410 G2M 3
## AAACCCATCCTAGCCT-1 Baf53Nb.IF2 -0.028925620 -0.0132582758 G1 0
## AAACCCATCGAAACAA-1 Baf53Nb.CL3 -0.008077289 -0.0028336129 G1 13
## seurat_clusters sort_clusters preAnno1 preAnno2 cnt
## AAACCCACAAGACGAC-1 3 3 PEMN.5 PEMN Baf53Nb.IF
## AAACCCAGTACAGGTG-1 0 0 PEMN.2 PEMN Baf53Nb.CL
## AAACCCAGTGGGCTCT-1 7 7 PSN.1 PSN Baf53Nb.CL
## AAACCCAGTTTGTTCT-1 3 3 PEMN.5 PEMN Baf53Nb.IF
## AAACCCATCCTAGCCT-1 0 0 PEMN.2 PEMN Baf53Nb.IF
## AAACCCATCGAAACAA-1 13 13 PIN.2 PIN Baf53Nb.CL
## pANN_0.25_0.005_471 DoubletFinder0.05 pANN_0.25_0.005_942
## AAACCCACAAGACGAC-1 0.07936508 Singlet 0.07936508
## AAACCCAGTACAGGTG-1 0.58730159 Doublet 0.58730159
## AAACCCAGTGGGCTCT-1 0.00000000 Singlet 0.00000000
## AAACCCAGTTTGTTCT-1 0.03174603 Singlet 0.03174603
## AAACCCATCCTAGCCT-1 0.01587302 Singlet 0.01587302
## AAACCCATCGAAACAA-1 0.00000000 Singlet 0.00000000
## DoubletFinder0.1
## AAACCCACAAGACGAC-1 Singlet
## AAACCCAGTACAGGTG-1 Doublet
## AAACCCAGTGGGCTCT-1 Singlet
## AAACCCAGTTTGTTCT-1 Singlet
## AAACCCATCCTAGCCT-1 Singlet
## AAACCCATCGAAACAA-1 Singlet
table(GEX.seur@meta.data[,c("preAnno1","DoubletFinder0.05")])
## DoubletFinder0.05
## preAnno1 Doublet Singlet
## PEMN.1 0 799
## PEMN.2 22 1076
## PEMN.3 22 697
## PEMN.4 27 458
## PEMN.5 5 739
## PEMN.6 16 287
## PIMN.1 4 923
## PIMN.2 6 567
## PIMN.3 6 561
## PIMN.4 8 223
## PIN.1 10 182
## PIN.2 25 243
## PIN.3 7 100
## PSN.1 36 457
## PSN.2 33 335
## PSN.3 20 204
## PSN.4 38 409
## PSN.5 7 40
## PSVN.1 1 39
## PSVN.2 11 146
## PSVN.3 26 256
## MIX.1 121 50
## MIX.2 14 99
## Glia 6 62
table(GEX.seur@meta.data[,c("preAnno1","DoubletFinder0.1")])
## DoubletFinder0.1
## preAnno1 Doublet Singlet
## PEMN.1 1 798
## PEMN.2 61 1037
## PEMN.3 56 663
## PEMN.4 64 421
## PEMN.5 13 731
## PEMN.6 44 259
## PIMN.1 14 913
## PIMN.2 24 549
## PIMN.3 27 540
## PIMN.4 36 195
## PIN.1 24 168
## PIN.2 47 221
## PIN.3 29 78
## PSN.1 52 441
## PSN.2 51 317
## PSN.3 31 193
## PSN.4 56 391
## PSN.5 20 27
## PSVN.1 16 24
## PSVN.2 25 132
## PSVN.3 38 244
## MIX.1 138 33
## MIX.2 42 71
## Glia 33 35
GEX.seur <- subset(GEX.seur, subset = DoubletFinder0.1=="Singlet" & preAnno2 %in% c("PEMN","PIMN","PIN","PSN","PSVN"))
GEX.seur
## An object of class Seurat
## 23258 features across 8342 samples within 1 assay
## Active assay: RNA (23258 features, 1500 variable features)
## 3 dimensional reductions calculated: pca, tsne, umap
DimPlot(GEX.seur, reduction = "umap", group.by = "preAnno1", label = T, label.size = 2.5) +
DimPlot(GEX.seur, reduction = "umap", group.by = "preAnno2", label = T, repel = F)
FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))
# Visualize QC metrics as a violin plot
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0.01, group.by = "FB.info", cols = color.FB)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0, group.by = "FB.info", cols = color.FB)
plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "FB.info", cols = color.FB)
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "FB.info", cols = color.FB)
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "FB.info", cols = color.FB)
plota + plotb + plotc
par(mar=c(6,4,2,3))
sl_stat <- table(GEX.seur$FB.info)
barplot(sl_stat,ylim = c(0,2300),
col = c("#FF6C91","lightgrey",color.FB),
main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:10*1.2-0.48 ,levels(GEX.seur$FB.info), las=3, cex.axis=0.85)
text(x=1:10*1.2-0.45,y=sl_stat+125,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)
#GEX.seur <- NormalizeData(GEX.seur, normalization.method = "LogNormalize", scale.factor = 10000)
GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1500)
# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2
head(VariableFeatures(GEX.seur), 100)
## [1] "Mgat4c" "Gal" "Zfp804a" "Grm5"
## [5] "Cntn5" "Nmu" "Klhl1" "Adgrg6"
## [9] "Gm32647" "4930432L08Rik" "Col24a1" "Cpne4"
## [13] "Kcnb2" "Gm30613" "Robo2" "Ntng1"
## [17] "Cntnap2" "Tmeff2" "Brinp3" "Gm29521"
## [21] "Nrxn3" "Lingo2" "Ano2" "Gpr149"
## [25] "Dgkg" "Cdh8" "Cdh18" "Zfp804b"
## [29] "Gm21847" "Ebf1" "Cdh9" "Ntrk3"
## [33] "Gm49953" "Gm15261" "Dgki" "Sgcz"
## [37] "Pdzrn4" "Cmah" "Nxph2" "Unc5d"
## [41] "Astn2" "Kcnip4" "1700111E14Rik" "Sema5a"
## [45] "Pcdh11x" "Schip1" "Pcdh9" "Prkg1"
## [49] "Itgb6" "Vwc2l" "Cdh6" "Nrg1"
## [53] "Postn" "Dcc" "Calcb" "4930509J09Rik"
## [57] "Car10" "Piezo2" "Grin3a" "Myl1"
## [61] "March1" "Clstn2" "Vip" "Gm38405"
## [65] "4930422I22Rik" "Rasgef1b" "Bmpr1b" "Efr3a"
## [69] "Ccbe1" "Rab27b" "Plxna4" "Fut9"
## [73] "Dpp10" "Gpc5" "Hbb-y" "AC150683.1"
## [77] "Gm49938" "Gm20754" "Gm30094" "Pbx3"
## [81] "Gm15680" "Asic2" "Egfem1" "9530059O14Rik"
## [85] "Pcdh10" "Abi3bp" "Cysltr2" "Gm15584"
## [89] "Serpini1" "Col25a1" "Skap1" "Nkain3"
## [93] "Robo1" "Fgf14" "Gm28494" "Trhde"
## [97] "Tafa2" "Spock3" "Chrna9" "Zfhx3"
GEX.seur <- ScaleData(GEX.seur, features = rownames(GEX.seur))
## Centering and scaling data matrix
# exclude MT genes and more
#DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^Gm|^Hist",rownames(GEX.seur),value = T)
DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^Gm|^Hist|Rik$|^AC|-ps",
rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
MT_gene <- grep("^mt-",rownames(GEX.seur),value = T)
GEX.seur <- RunPCA(GEX.seur, features = setdiff(VariableFeatures(object = GEX.seur),
c(MT_gene,
DIG,
CC_gene) ))
## PC_ 1
## Positive: Ntrk3, Tmeff2, Robo2, Cdh8, Ano2, Nrxn3, Cpne4, Clstn2, Myl1, Mgat4c
## Dgkg, Gpr149, Pdzrn4, Adgrg6, Plxna4, Grin3a, Pcdh10, Cdh6, Spock3, Itgb6
## Cysltr2, Cntn5, Zfp804a, Sgcz, Iqgap2, Cacna2d3, Cntnap2, Ccbe1, Astn2, Cux2
## Negative: Lrrtm4, Galntl6, Tox, Ndst4, Cacnb2, Grik1, Kcnd2, Lrp1b, Bnc2, Pcdh7
## Epha5, Synpo2, Kcnab1, Syt6, Tshz2, Zfp536, Plcxd3, Cdc14a, Man2a1, Chrna7
## Specc1, Lrrc4c, Nos1, Galnt17, Pitpnc1, Lrrc7, Tenm3, Fat3, Dach1, Brinp2
## PC_ 2
## Positive: Rbfox1, Bnc2, Ptprt, Tafa1, Tshz2, Grik1, Gpc6, Frmd4b, Plcxd3, St6galnac3
## Tox, Caln1, Brinp2, Fbxw15, Dock2, Cdc14a, Tcf7l2, Fbxw24, Tmem132c, Pcdh7
## Pld5, Oprk1, Sdk2, Specc1, Tafa5, Adamtsl1, Slc35f4, Nfia, Unc5c, Stxbp5l
## Negative: Nos1, Auts2, Etv1, Egfem1, Cadps2, Kcnq5, Gfra1, Schip1, Asic2, Dgkb
## Cmah, Dach1, Kcnt2, Epha5, Kcnab1, Rgs6, Stxbp6, Alcam, Creb5, Cntnap5a
## Il1rapl1, Lrrc4c, Hs6st3, Tmem108, Adgrl3, Cdh11, Gria3, Ebf1, Gramd1b, Adcy2
## PC_ 3
## Positive: Cdh18, Klhl1, Pbx3, Kcnip4, Meis1, Sema5a, Kctd8, C79798, Gabrg3, Htr4
## Vwc2l, Gpc6, Cntn3, Dlc1, Serpini1, Zbbx, Csmd3, Zfhx3, Cdh9, March1
## Galnt18, Skap1, Plcl1, Stxbp5l, Pakap, Pde4d, Khdrbs2, Csmd1, Cadm2, Nrp2
## Negative: Adgrg6, Sgcd, Cysltr2, Slc4a4, Itgb6, Nfia, Nmu, Nos1, Gpr149, Grin3a
## Dapk2, Ano2, Rgs6, Cbln2, Lhfpl2, Zfp804a, Ccbe1, Dgkg, Zfp536, Kcnab1
## Otof, Efr3a, Cntnap5a, Cpne4, Zeb2, Smad6, Syt15, Gfra1, Scn11a, Zbtb7c
## PC_ 4
## Positive: Dock10, Kcnip4, Lingo2, Csmd3, Sorcs1, Gda, Ndst4, Cntn5, Tac1, Nxph1
## Cntn3, Thsd4, Kctd8, Nyap2, Sgcz, Kcnq5, Lrp1b, Robo1, Plcb1, Pld5
## Tenm2, Lin7a, Lrrc7, Abca5, Lrrtm4, Dmd, Pgm5, Kctd16, Sorcs3, Prkg1
## Negative: Klhl1, Sema5a, Vwc2l, March1, Rasgef1b, Cdh9, Alk, Prune2, Galnt18, Il1rapl1
## Chsy3, Mgat4c, Kcnh7, Pbx3, Bcl2, Zfhx3, Dpp10, Zbbx, Galnt17, Grm5
## Pcdh7, C79798, Vcan, Kcnb2, Lncbate1, Auts2, Ghr, Dcc, Scgn, Olfr78
## PC_ 5
## Positive: Prkg1, Kcnt2, Dgkb, Alcam, Car10, Rasgef1b, Galntl6, Epha5, Mgat4c, Vwc2l
## Vcan, Thsd7b, Cdh20, Dmd, Dpp10, Rgs6, Sema5a, Gucy1a2, Lingo2, Gabrg3
## Cdh9, Khdrbs2, Slc2a13, Kctd8, Lrrc4c, Htr4, Antxr2, Ptprz1, Adcy2, Man1a
## Negative: Ntng1, Ebf1, Trps1, Trhde, Cntn4, Csmd1, Gna14, Zmat4, Col18a1, Nkain3
## Tmtc2, Myo1e, Tenm4, Sctr, Nxph2, Cpa6, Nrg1, Kcnd2, Npas3, Neurod6
## Fstl4, Sez6l, Nav2, Ccdc85a, Myo16, Gal, Glp1r, Fbn2, Tfcp2l1, Nrp1
length(setdiff(VariableFeatures(object = GEX.seur),
c(MT_gene,DIG, CC_gene) ))
## [1] 1067
head(setdiff(VariableFeatures(object = GEX.seur),
c(MT_gene,DIG, CC_gene) ),300)
## [1] "Mgat4c" "Gal" "Zfp804a" "Grm5" "Cntn5"
## [6] "Nmu" "Klhl1" "Adgrg6" "Col24a1" "Cpne4"
## [11] "Kcnb2" "Robo2" "Ntng1" "Cntnap2" "Tmeff2"
## [16] "Brinp3" "Nrxn3" "Lingo2" "Ano2" "Gpr149"
## [21] "Dgkg" "Cdh8" "Cdh18" "Zfp804b" "Ebf1"
## [26] "Cdh9" "Ntrk3" "Dgki" "Sgcz" "Pdzrn4"
## [31] "Cmah" "Nxph2" "Unc5d" "Astn2" "Kcnip4"
## [36] "Sema5a" "Pcdh11x" "Schip1" "Pcdh9" "Prkg1"
## [41] "Itgb6" "Vwc2l" "Cdh6" "Nrg1" "Postn"
## [46] "Dcc" "Calcb" "Car10" "Piezo2" "Grin3a"
## [51] "Myl1" "March1" "Clstn2" "Vip" "Rasgef1b"
## [56] "Bmpr1b" "Efr3a" "Ccbe1" "Rab27b" "Plxna4"
## [61] "Fut9" "Dpp10" "Gpc5" "Pbx3" "Asic2"
## [66] "Egfem1" "Pcdh10" "Abi3bp" "Cysltr2" "Serpini1"
## [71] "Col25a1" "Skap1" "Nkain3" "Robo1" "Fgf14"
## [76] "Trhde" "Tafa2" "Spock3" "Chrna9" "Zfhx3"
## [81] "Csmd3" "Galnt13" "Ghr" "Tac1" "Pde4d"
## [86] "Luzp2" "Pde1a" "Tbx20" "Trps1" "Kcnq5"
## [91] "Mid1" "Cemip" "Dgkb" "Myo3b" "Kcnh7"
## [96] "Acp7" "Nos1" "Flrt2" "Plcl1" "Iqgap2"
## [101] "Cacna2d3" "Csmd1" "Tnr" "Kctd16" "Cadm2"
## [106] "Parm1" "Gna14" "Cntn4" "Chrm3" "Nell1"
## [111] "Abca9" "Slc4a4" "Il1rapl1" "Nxph1" "Kctd8"
## [116] "Adamts6" "Fbxl7" "Myh3" "Cartpt" "Hs6st3"
## [121] "Arhgap6" "Kif26b" "Penk" "Aff2" "Grm7"
## [126] "Galnt18" "Cacnb4" "Synpr" "Capn9" "Sst"
## [131] "Lama2" "Ano5" "Cdh10" "Grik1" "Opcml"
## [136] "Kcnd2" "Cadps2" "Exoc3l2" "Esr1" "Galntl6"
## [141] "L3mbtl4" "Prune2" "Gabrg3" "Spag16" "Epha7"
## [146] "Cntn3" "Prr16" "Cux2" "Calca" "Chsy3"
## [151] "Dach1" "Antxr2" "Morc1" "Met" "Thsd7b"
## [156] "Kirrel3" "Hs3st4" "Zbbx" "Npas3" "Terb2"
## [161] "Cbln2" "Loxl1" "Il18r1" "Stab1" "Galr1"
## [166] "Runx1" "Egfr" "Ntrk2" "Dock10" "Olfr78"
## [171] "Rasgrf1" "Scnn1a" "Scn7a" "Cpa6" "Col8a1"
## [176] "Shisa9" "D5Ertd615e" "C79798" "Pde7b" "Rerg"
## [181] "Terf1" "Ttc29" "Pak7" "Vcan" "Etv1"
## [186] "Grp" "Olfr889" "Adam2" "Rarb" "Tenm2"
## [191] "Lhfpl2" "Sorcs3" "Gabrb1" "Alcam" "Bpifc"
## [196] "Otof" "Necab1" "Gng2" "Calcrl" "Rxfp2"
## [201] "Dapk2" "Slco1a1" "Ntm" "Pdgfra" "Plxdc2"
## [206] "Lrrc4c" "Fstl4" "Bche" "Mrc1" "Nell1os"
## [211] "Stxbp6" "Hcn1" "F13a1" "Htr4" "Nr4a3"
## [216] "Elf5" "Slc7a15" "Myo1h" "Fzd7" "Thsd4"
## [221] "Creb5" "Sorcs1" "Lockd" "Tenm4" "Tafa1"
## [226] "Dysf" "Rxfp1" "Nwd2" "Igfbp5" "Rmst"
## [231] "Gfra1" "Dlc1" "Ppp3ca" "Plcb1" "Rtl4"
## [236] "Cfh" "Cblc" "Sctr" "Arhgap15" "Nfatc1"
## [241] "Rgs6" "Bmp5" "Slc44a5" "Ngfr" "Ifi203"
## [246] "Lamc3" "Itgbl1" "Nox4" "Cntnap5b" "Naaladl2"
## [251] "Tex15" "Htr2c" "Adgrl2" "Scube1" "Rbpms"
## [256] "Tex21" "Casp8" "Syt9" "Etl4" "Dab2"
## [261] "Ptprt" "Snx7" "Tll1" "Moxd1" "Plac8"
## [266] "Grid1" "Agrn" "Prkag2" "Bcl2" "Tenm1"
## [271] "Prkca" "C1qtnf7" "Ror1" "Dab1" "Eya1"
## [276] "Pde1c" "Edil3" "Hs3st2" "Fam47e" "Fbxw21"
## [281] "Dgat2l6" "Col18a1" "Tmem64" "Trpc3" "Adam12"
## [286] "Cfap299" "Best2" "Hsd17b14" "Cpne8" "Vit"
## [291] "Adcy2" "Myo1e" "Cldn5" "Rad51b" "Brinp2"
## [296] "Fras1" "Fgf10" "Zbtb7c" "Grm8" "Atpaf2"
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB) +
DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB)
DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)
ElbowPlot(GEX.seur,ndims = 50)
PCs <- 1:18
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param =20)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 1.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8342
## Number of edges: 333410
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8336
## Number of communities: 21
## Elapsed time: 0 seconds
GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity=100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 133)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 15:21:35 UMAP embedding parameters a = 0.9922 b = 1.112
## 15:21:35 Read 8342 rows and found 18 numeric columns
## 15:21:35 Using Annoy for neighbor search, n_neighbors = 20
## 15:21:35 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 15:21:36 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\Rtmp6hZbq2\file787045932716
## 15:21:36 Searching Annoy index using 1 thread, search_k = 2000
## 15:21:37 Annoy recall = 100%
## 15:21:38 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 20
## 15:21:39 Initializing from normalized Laplacian + noise (using irlba)
## 15:21:39 Commencing optimization for 500 epochs, with 239650 positive edges
## 15:22:00 Optimization finished
DimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)
DimPlot(GEX.seur, reduction = "umap", group.by = "preAnno1", label = T, label.size = 2.5) +
DimPlot(GEX.seur, reduction = "umap", group.by = "preAnno2", label = T, repel = F)
FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))
DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "FB.info", split.by = "FB.info",
ncol = 4, cols = color.FB)
DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "FB.info", split.by = "FB.info",
ncol = 4, cols = color.FB)
ref.seur <- readRDS("../analysis_ref/GSE149524.P21.integration_Anno.s.rds")
ref.seur
## An object of class Seurat
## 37583 features across 4419 samples within 3 assays
## Active assay: SCT (16365 features, 0 variable features)
## 2 other assays present: RNA, integrated
## 3 dimensional reductions calculated: pca, tsne, umap
ref.seur <- RunUMAP(ref.seur, assay = "integrated", umap.method = "uwot-learn", dims = 1:30)
## Warning in RunUMAP.default(object = data.use, reduction.model =
## reduction.model, : 'uwot-learn' is deprecated. Set umap.method = 'uwot' and
## return.model = TRUE
## UMAP will return its model
## 15:22:15 UMAP embedding parameters a = 0.9922 b = 1.112
## 15:22:15 Read 4419 rows and found 30 numeric columns
## 15:22:15 Using Annoy for neighbor search, n_neighbors = 30
## 15:22:15 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 15:22:15 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\Rtmp6hZbq2\file7870317045c9
## 15:22:15 Searching Annoy index using 1 thread, search_k = 3000
## 15:22:16 Annoy recall = 100%
## 15:22:17 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 15:22:17 Initializing from normalized Laplacian + noise (using irlba)
## 15:22:17 Commencing optimization for 500 epochs, with 181154 positive edges
## 15:22:30 Optimization finished
anchors.ref <- FindTransferAnchors(
reference = ref.seur,
reference.assay = "integrated",
query = GEX.seur,
#query.assay = "RNA",
query.assay = "RNA",
reference.reduction = 'pca',
normalization.method = "LogNormalize",
#normalization.method = "SCT",
reduction = "pcaproject",
dims = 1:30,
k.anchor = 80,
k.filter = 120,
k.score = 60,
max.features = 1600,
nn.method = "annoy",
n.trees = 100,
eps = 0
)
## Projecting cell embeddings
## Finding neighborhoods
## Finding anchors
## Found 20720 anchors
## Filtering anchors
## Retained 6498 anchors
local_neur.pred <- MapQuery(
anchorset = anchors.ref,
query = GEX.seur,
reference = ref.seur,
refdata = list(
ref.Anno1 = 'Anno1',
ref.Anno2 = 'Anno2'
),
reference.reduction = 'pca',
reduction.model = 'umap'
)
## Warning: `invoke()` is deprecated as of rlang 0.4.0.
## Please use `exec()` or `inject()` instead.
## This warning is displayed once per session.
## Finding integration vectors
## Finding integration vector weights
## Predicting cell labels
## Warning: Keys should be one or more alphanumeric characters followed
## by an underscore, setting key from predictionscoreref.Anno1_ to
## predictionscorerefAnno1_
## Predicting cell labels
## Warning: Keys should be one or more alphanumeric characters followed
## by an underscore, setting key from predictionscoreref.Anno2_ to
## predictionscorerefAnno2_
##
## Integrating dataset 2 with reference dataset
## Finding integration vectors
## Integrating data
## Computing nearest neighbors
## Running UMAP projection
## 15:24:14 Read 8342 rows
## 15:24:14 Processing block 1 of 1
## 15:24:14 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 15:24:14 Initializing by weighted average of neighbor coordinates using 1 thread
## 15:24:14 Commencing optimization for 167 epochs, with 250260 positive edges
## 15:24:22 Finished
local_neur.pred
## An object of class Seurat
## 23284 features across 8342 samples within 3 assays
## Active assay: RNA (23258 features, 1500 variable features)
## 2 other assays present: prediction.score.ref.Anno1, prediction.score.ref.Anno2
## 5 dimensional reductions calculated: pca, tsne, umap, ref.pca, ref.umap
#
raw1_p1 <- DimPlot(local_neur.pred,
reduction = 'umap',
group.by = 'preAnno1',
label=TRUE,
label.size = 3,
repel = F) + NoLegend()+ labs(title="Baf53Nb preAnno1")
raw1_p2 <- DimPlot(local_neur.pred,
reduction = 'umap',
group.by = 'predicted.ref.Anno1',
label=TRUE,
label.size = 3,
repel = T) + NoLegend()
raw1_p3 <- DimPlot(local_neur.pred,
reduction = 'umap',
group.by = 'preAnno2',
label=TRUE,
label.size = 3,
repel = TRUE) + NoLegend()+ labs(title="Baf53Nb preAnno2")
raw1_p4 <- DimPlot(local_neur.pred,
reduction = 'umap',
group.by = 'predicted.ref.Anno2',
label=TRUE,
label.size = 3,
repel = TRUE) + NoLegend()
cowplot::plot_grid(raw1_p1,raw1_p2,ncol = 2)
cowplot::plot_grid(raw1_p3,raw1_p4,ncol = 2)
cowplot::plot_grid(FeaturePlot(local_neur.pred, features = "predicted.ref.Anno1.score", reduction = "umap"),
FeaturePlot(local_neur.pred, features = "predicted.ref.Anno2.score", reduction = "umap"),ncol = 2)
clusters sort new orders for mapping comparison,
and do newAnno based on it through ref.mapping to NatNeur2021,
EMN1-5 C0/1/3-C9-C2-C14-C16
IMN1-3 C5/4/15-C6-C11
IN1-3 C12-C18-C20(Sst+)
IPAN1-4 C7/10-C8/17-C19-C13
GEX.seur$sort_new <- factor(as.character(GEX.seur$seurat_clusters),
levels = c(0,1,3,
9,
2,
14,
16,
5,4,15,
6,
11,
12,
18,
20,
7,10,
8,17,
19,
13))
#
local_neur.pred$sort_new <- GEX.seur$sort_new
pred.ref.Anno1 <- table(prediction=local_neur.pred$predicted.ref.Anno1,
clusters=local_neur.pred$sort_new)[paste0("ENC",c(1,#2,
3,4,
8,9,
10,#11,5,
6,7,12)),]
pred.ref.Anno1
## clusters
## prediction 0 1 3 9 2 14 16 5 4 15 6 11 12 18 20 7 10
## ENC1 943 884 683 288 125 12 4 63 120 0 185 145 0 1 24 1 0
## ENC3 0 0 0 106 586 151 0 0 0 0 0 0 0 0 0 0 0
## ENC4 0 0 0 0 0 55 200 0 0 0 0 0 0 0 0 0 0
## ENC8 0 0 0 0 0 0 0 470 449 213 227 0 0 1 0 0 0
## ENC9 0 0 0 0 0 0 0 0 1 0 63 123 0 0 0 0 0
## ENC10 0 1 0 0 28 13 2 15 90 0 10 20 257 134 0 0 0
## ENC6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 423 316
## ENC7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 11
## ENC12 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0
## clusters
## prediction 8 17 19 13
## ENC1 0 1 9 2
## ENC3 0 0 0 0
## ENC4 0 0 0 0
## ENC8 0 0 1 0
## ENC9 0 0 0 0
## ENC10 0 0 7 3
## ENC6 0 0 0 0
## ENC7 405 178 1 0
## ENC12 0 0 45 230
pred.ref.Anno2 <- table(prediction=local_neur.pred$predicted.ref.Anno2,
clusters=local_neur.pred$sort_new)
pred.ref.Anno2
## clusters
## prediction 0 1 3 9 2 14 16 5 4 15 6 11 12 18 20 7 10
## EMN1 943 884 683 290 197 25 4 63 120 0 185 145 0 1 24 1 0
## EMN3 0 0 0 104 287 10 0 0 0 0 0 0 0 0 0 0 0
## EMN4 0 0 0 0 222 111 0 0 0 0 0 0 0 0 0 0 0
## EMN5 0 0 0 0 0 72 200 0 0 0 0 0 0 0 0 0 0
## IMN1 0 0 0 0 0 0 0 470 449 213 227 0 0 1 0 0 0
## IMN2 0 0 0 0 0 0 0 0 1 0 63 123 0 0 0 0 0
## IN1 0 1 0 0 33 13 2 15 90 0 10 20 257 134 0 0 0
## IN2 0 0 0 0 0 0 0 0 0 0 0 0 0 7 0 0 0
## IPAN1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 423 316
## IPAN2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 11
## IPAN4 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 0 0
## clusters
## prediction 8 17 19 13
## EMN1 0 1 9 2
## EMN3 0 0 0 0
## EMN4 0 0 0 0
## EMN5 0 0 0 0
## IMN1 0 0 1 0
## IMN2 0 0 0 0
## IN1 0 0 7 3
## IN2 0 0 0 0
## IPAN1 0 0 0 0
## IPAN2 405 178 1 0
## IPAN4 0 0 45 230
cowplot::plot_grid(
pheatmap::pheatmap(pred.ref.Anno1,
main = "Cell Count (Baf53Nb mapping to NatNeur2021 P21, Anno1)",
gaps_row = c(3,5,6),
gaps_col = c(7,12,15),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(t(100*rowRatio(t(pred.ref.Anno1))),
main = "Cell Ratio (Baf53Nb mapping to NatNeur2021 P21, Anno1)",
gaps_row = c(3,5,6),
gaps_col = c(7,12,15),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent = T)$gtable,
ncol = 1)
cowplot::plot_grid(
pheatmap::pheatmap(pred.ref.Anno2,
main = "Cell Count (Baf53Nb mapping to NatNeur2021 P21, Anno2)",
gaps_row = c(4,6,8),
gaps_col = c(7,12,15),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(t(100*rowRatio(t(pred.ref.Anno2))),
main = "Cell Ratio (Baf53Nb mapping to NatNeur2021 P21, Anno2)",
gaps_row = c(4,6,8),
gaps_col = c(7,12,15),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent = T)$gtable,
ncol = 1)
GEX.seur$newAnno <- as.character(GEX.seur$seurat_clusters)
GEX.seur$newAnno[GEX.seur$newAnno %in% c(0,1,3)] <- "EMN1"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(9)] <- "EMN2"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(2)] <- "EMN3"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(14)] <- "EMN4"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(16)] <- "EMN5"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(5,4,15)] <- "IMN1"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(6)] <- "IMN2"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(11)] <- "IMN3"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(12)] <- "IN1"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(18)] <- "IN2"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(20)] <- "IN3"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(7,10)] <- "IPAN1"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(8,17)] <- "IPAN2"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(19)] <- "IPAN3"
GEX.seur$newAnno[GEX.seur$newAnno %in% c(13)] <- "IPAN4"
GEX.seur$newAnno <- factor(GEX.seur$newAnno,
levels = c(paste0("EMN",1:5),
paste0("IMN",1:3),
paste0("IN",1:3),
paste0("IPAN",1:4)))
# define newAnno colors
color.newA <- c("#8AB6CE","#678BB1","#3975C1","#669FDF","#4CC1BD",
"#BF7E6B","#D46B35","#FF8080",
"#BDAE8D","#BD66C4","#C03778",
"#97BA59","#DFAB16","#2BA956","#9FE727")
cowplot::plot_grid(
DimPlot(GEX.seur, reduction = "umap", group.by = "seurat_clusters", label = T, label.size = 3.25,repel = F, pt.size = 0.45),
DimPlot(GEX.seur, label = F, pt.size = 0.45, repel = F, reduction = 'umap', group.by = "cnt",
ncol = 1, cols = color.FB[c(5,1)]) ,
rel_widths = c(4.8,5),
ncol = 2)
cowplot::plot_grid(
DimPlot(GEX.seur, reduction = "umap", group.by = "newAnno", label = T, label.size = 3.25,repel = F, pt.size = 0.45),
DimPlot(GEX.seur, label = F, pt.size = 0.45, repel = F, reduction = 'umap', group.by = "cnt",
ncol = 1, cols = color.FB[c(5,1)]) ,
rel_widths = c(4.75,5),
ncol = 2)
cowplot::plot_grid(
DimPlot(GEX.seur, reduction = "umap", group.by = "newAnno", label = T, label.size = 3.25,repel = F, pt.size = 0.45,
cols = color.newA),
DimPlot(GEX.seur, label = F, pt.size = 0.45, repel = F, reduction = 'umap', group.by = "cnt",
ncol = 1, cols = color.FB[c(5,1)]) ,
rel_widths = c(4.75,5),
ncol = 2)
GEX.seur@meta.data$cnt <- factor(gsub("1$|2$|3$|4$","",as.character(GEX.seur@meta.data$FB.info)))
GEX.seur$rep <- paste0("rep",
gsub("Baf53Nb.IF|Baf53Nb.CL","",as.character(GEX.seur$FB.info)))
stat_newAnno <- GEX.seur@meta.data[,c("cnt","rep","newAnno")]
stat_newAnno.s <- stat_newAnno %>%
group_by(cnt,rep,newAnno) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup
## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.
pcom.newAnno <- stat_newAnno.s %>%
ggplot(aes(x = newAnno, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = c("skyblue","pink"), name="") +
labs(title="Cell Composition of Baf53Nb data, newAnno", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=newAnno, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom.newAnno
glm.nb - anova.LRT
Sort.N <- c("Baf53Nb.CL","Baf53Nb.IF")
glm.nb_anova.newAnno <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat_newAnno.s_N <- stat_newAnno.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat_newAnno.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat_newAnno.s_N$total <- total.N[paste0(stat_newAnno.s_N$cnt,
"_",
stat_newAnno.s_N$rep),"total"]
glm.nb_anova.newAnno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat_newAnno.s_N$newAnno)){
glm.nb_anova.newAnno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat_newAnno.s_N %>% filter(newAnno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat_newAnno.s_N %>% filter(newAnno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat_newAnno.s_N %>% filter(newAnno==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.newAnno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.newAnno[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}
glm.nb_anova.newAnno_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.newAnno)))
rownames(glm.nb_anova.newAnno_df) <- names(glm.nb_anova.newAnno)
colnames(glm.nb_anova.newAnno_df) <- gsub("X","C",colnames(glm.nb_anova.newAnno_df))
glm.nb_anova.newAnno_df
## EMN1 EMN2 EMN3 EMN4 EMN5
## Baf53Nb.CLvsBaf53Nb.IF 0.008660417 0.3309766 0.891037 0.8971586 0.6479739
## IMN1 IMN2 IMN3 IN1 IN2
## Baf53Nb.CLvsBaf53Nb.IF 0.07575576 0.1147984 0.09201819 0.8945842 0.5844769
## IN3 IPAN1 IPAN2 IPAN3 IPAN4
## Baf53Nb.CLvsBaf53Nb.IF 0.8401056 0.09909844 0.4667117 0.9061248 0.7642453
markers.new <- list(EMN=c("Chat","Gfra2","Casz1","Xylt1",
"Ptprt","Bnc2","Tox","Oprk1",
"Kalrn","Pi15","Drd2","Adamtsl1",
"Fbxw15","Fbxw24","Cdc14a","Chrna7",
"Caln1","Tmem132c","Satb1","Cntnap5b",
"Gabrb1","Nxph1","Lama2","Lrrc7",
"Ryr3","Eda","Chgb","Pgm5",
"Shc4","Vgll3","Ptn","Tac1",
"Kctd8","Ntrk2","Penk","Pde7b",
"Fut9","Nfatc1","Egfr","Mgll",
"Ntn1","Prlr","Chrm3"
),
IMN=c("Nos1","Dach1","Kcnab1",
"Gfra1","Etv1",
"Man1a","Airn","Enpp1","Unc13c",
"Plpp3","Fat1","Adcy2","S100a4",
"Npy","Cmah","Creb5","Vip","Pde1a",
"Ebf1","Gpc5","Col25a1","Cartpt",
"Lgr5","Gabrb2"
),
IN=c("Npas3","Synpr","St18","Gal",
"Kcnk13",
"Neurod6","Moxd1","Sctr","Fndc1",
"Piezo1","Sst","Adamts9","Kcnn2",
"Pantr1","Vwc2","Vipr2","Tacr1",
"Calb2"),
IPAN=c("Calcb","Nmu","Kcnj12","Nog",
"Bmp4","Adgrg6","Cysltr2","Pcdh10",
"Ngfr","Galr1","Met",
"Htr3a","Il7","Aff2","Gpr149",
"Efr3a","Cdh6","Cdh8","Pdzrn4",
"Clstn2","Cachd1","Ano2","Ntrk3",
"Cpne4","Vwc2l","Cdh9","Scgn",
"Vcan","Cck","Kcnh7","Piezo2",
"Abca6","Fam107b","Npy1r","Abca9",
"Abca8b","Rerg","Bmpr1b","Skap1",
"L3mbtl4","Tafa2","Nxph2","Gm32647",
"Gm29521","Ntng1"))
pm.CL_new <- DotPlot(subset(GEX.seur,subset=cnt %in% c("Baf53Nb.CL")), features = as.vector(unlist(markers.new)), group.by = "newAnno",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CL only")
pm.CL_new
pm.All_new <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new)), group.by = "newAnno",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new
markers.new.s <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
"Gfra2","Oprk1","Adamtsl1",
"Fbxw15","Fbxw24","Chrna7",
"Satb1","Cntnap5b",
"Gabrb1","Nxph1","Lama2","Lrrc7",
"Ryr3","Eda","Tac1",
"Kctd8","Ntrk2","Penk",
"Fut9","Nfatc1","Egfr","Mgll",
"Chrm3"
),
IMN=c("Nos1","Kcnab1",
"Gfra1","Etv1",
"Man1a","Airn",
"Adcy2","Cmah","Creb5","Vip","Pde1a",
"Ebf1","Gpc5"
),
IN=c("Npas3","Synpr","St18","Gal",
"Kcnk13",
"Neurod6","Moxd1","Sctr",
"Piezo1","Sst","Adamts9","Kcnn2",
"Calb2"),
IPAN=c("Calcb","Nmu","Adgrg6","Pcdh10",
"Ngfr","Galr1","Il7","Aff2",
"Gpr149","Cdh6","Cdh8",
"Clstn2","Ano2","Ntrk3",
"Cpne4","Vwc2l","Cdh9","Scgn",
"Vcan","Cck","Piezo2","Kcnh7",
"Rerg","Bmpr1b","Skap1","Ntng1",
"Tafa2","Nxph2"))
pm.CL_new.s <- DotPlot(subset(GEX.seur,subset=cnt %in% c("Baf53Nb.CL")), features = as.vector(unlist(markers.new.s)), group.by = "newAnno",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="CL only")
pm.CL_new.s
pm.All_new.s <- DotPlot(GEX.seur, features = as.vector(unlist(markers.new.s)), group.by = "newAnno",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All")
pm.All_new.s
#saveRDS(GEX.seur,"./Baf53Nb.Ileum_IFd5.pure_newAnno20230728.rds")